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ABSTRACT 

In a series of two papers, we test the primordial scenario of globular cluster formation using results of high- 
resolutions iV-body simulations. In this first paper we study the initial relaxation of a stellar core inside a live 
dark matter minihalo in the early universe. Our dark-matter dominated globular clusters show features which 
are usually attributed to the action of the tidal field of the host galaxy. Among them are the presence of an 
apparent cutoff ("tidal radius") or of a "break" in the outer parts of the radial surface brightness profile, and a 
flat line-of-sight velocity dispersion profile in the outskirts of the cluster. The apparent mass-to-light ratios of 
our hybrid (stars + dark matter) globular clusters are very close to those of purely stellar clusters. We suggest 
that additional observational evidence such as the presence of obvious tidal tails is required to rule out the 
presence of significant amounts of dark matter in present day globular clusters. 

Subject headings: globular clusters: general — methods: A^-body simulations — dark matter — early universe 



1. INTRODUCTION 

The origin of globular clusters (GCs) is one of the important 
unsolved astrophysical problems. One can divide all proposed 
scenarios for GC formation into two large categories: 1) in 
situ formation (when GCs are formed inside their present day 
host galaxies), and 2) pregalactic formation (when GCs are 
formed in smaller galaxies which later merge to become a part 
of the present day host galaxy). 

Pregalactic scenarios for GC formation can explain many 
puzzling properties of GCs in the Milky Way and other galax- 
ies. This is especially true for so-called metal-poor GCs 
(with [Fe/H] < — 1), which form a distinctive population in 
many galaxies. Many of their properties are highly suggestive 
of their pregalactic origin: low metallicity, large age com- 
parable to the age of the universe, isotropic orbits with the 
apocentric-to-pericentric distances rati o of R a /R p ~ 5 .4 ± 3 .7 
JPinescu. Girard. & van Altenalll999l based on data for 38 
metal-poor Galactic GCs) which is consistent with the orbits 
of dark matter (DM) subhalos (R a /R p ~ 5 ... 6) in cosmolog- 
ical ACDM simulations, and a very weak or no correlation 
with the properties of the host galaxy. 

One very interesting variation of the pregalactic picture is 
the primordial scena rio of GC form ation, which was first pro- 
posed by Peebles & Dicke ( 1968). In its modern interpreta- 
tion, the primordial scenario assumes that GCs formed at the 
center of small DM minihalos in the early universe, before or 
during the reionization of the universe, which is believed to 
have ended around the redshift of z = 6 tBecker et all 2001 ). 
The primordial pict ure ha s been con sidered b y many a uthors, 
including Peebles (1984), Rosenblatt. Faber. & Blumenthal 
(119881). iPadoan Jimenez. & Jones! d!997l) ICenl (120011) . 
Bromm & Clarke! (120021) . and lBeaslev et all (120031) . 

It is usually assumed that during the reionization of the 
universe no star formation can take place inside small DM 
halos with the virial temperature below ~ 10 4 K (or virial 
mass less than ~ 10 s M Q ), as the gas in such halos should 
escape the shallow gravitational potential after being heated 
to ~ 10 4 K by the metagalactic Ly man continuum ( LyC) 
background (Barkana & Loeblll999t). Some authors (ICenl 
l200U lRicotti Gnedin. & Shul]ll2002l) suggested though, that 
the reionization of the universe can actually trigger star forma- 
tion in small DM minihalos through different positive feed- 



back mechanisms. In the scenario of Cen (2001), star forma- 
tion in small gas-rich minihalos is triggered by radially con- 
verging radiation shock fronts caused by the external LyC ra- 
diation fie ld. In the cosmolog ical simulations with radiative 
transfer of Ric otti et alJ (120021) . small-halo objects constitute 
the bulk of mass in stars until at least redshift z ~ 10 because 
of the increased non-equilibrium fraction of free electrons in 
front of the cosmological H II regions and inside the relic cos- 
mological H II regions, which results in more efficient H2 for- 
mation and hence cooling of the gas. 

Despite the fact that many authors considered different hy- 
drodynamic and radiative processes which can lead to forma- 
tion of GC-like stellar clusters inside DM minihalos, there 
has been no detailed study on what happens to such hybrid 
(GC + DM halo) object after the formation of stars from 
the point of galactic dynamics. Fully consistent cosmologi- 
cal sim ulations of str ucture formation in the universe (such as 
those of Ricotti et al. 2002j), and even higher resolution "semi- 
consistent" simulations (like in Bromm & Clarke 2002) lack 
orders of magnitude in spatial and mass resolution to be able 
to answer the following questions: 1) How does the presence 
of DM modify the observable properties of GCs? 2) Are there 
observable features from which the presence of DM in a GC 
can be inferred? 3) Will DM in "hybrid GCs" survive tidal 
stripping during the hierarchical assemblage of substructure 
leading to the formation of large galaxies such as Milky Way? 

We try to answer the above questions in a series of two pa- 
pers. In this first paper, we address the first and the second 
questions, with the last ques tion bei ng dealt with in the second 
paper ( Mashchenko & Sills 2004b [ hereafter Paper II). Using 
the N-body tree-code GADGET (Springel. Yoshida &Whitd 
12001 . we follow the relaxation of initially non-equilibrium 
stellar clusters inside live DM minihalos around the red- 
shift of z - 7. Our simulations are collisionless, and can 
be directly compared with dynamically unevolved GCs; 
the impact of secular evolution (core collapse) on our re- 
sults will be explored in Pape r II. DM halos have either 
iNavarro. Frenk. & WMtel i 19971 hereafter NFW) or iBurkertl 
7l995l) profiles. and have structural parameters taken from 
cosmological simulations. We use proto-GC model of 
Mashchenko & Sills (12004a! hereafter MS04) to set up the 
initial non-equilibrium configuration of stellar clusters inside 
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DM halos. In MS04 we showed that the collapse of homo- 
geneous isothermal stellar spheres leads to the formation of 
clusters with the surface brightness profiles being very sim- 
ilar to those of dynamically young Galactic GCs. In this 
model, all the observed correlations between structural and 
dynamic parameters of Galactic GCs are accurately repro- 
duced if the initial stellar density p,> and velocity dispersion 



er,.» had the same universal values: pi .* ~ 14 M Q pc and 
cr (> ~ 1.91 kms" 1 . 

This paper is organized as follows. In Section[2]we describe 
our method of simulating the evolution of a "hybrid" GC and 
list the physical and numerical parameters of our models. In 
Section|5]we show the results of the simulations. Finally, in 
Section|4]we discuss the results and give our conclusions. 

2. MODEL 

2.1. Initial Considerations 

In our model, GCs with a stellar mass m* are formed at 
the center of DM halos with the virial mass wdm at the red- 
shift of z = 7. We fix the mass ratio \ = 'm*/'«dm to 0.0088. 
The adopted value of x is between the unive rsal bary onic-to- 
DM density ratio U b /^dm = 0.20 ( Spe reel et alJ200 3 ) and the 
fraction of baryons in GCs in the modern universe of ~ 0.0025 
jMcLauehlinl 199 9). This leaves enough room for such effects 
as a less than 100% efficiency of star formation in proto-GCs, 
mass losses due to stellar evolution (through supernovae and 
stellar winds), a decay of GC systems due to dynamic evolu- 
tion of the clusters in the presence of tidal fields, and any bi- 
ased mechanism of GC formation (e.g., when GCs are formed 
only in DM halos with a low specific angular momentum, like 
in lCenl2001l) . 

2.2. DM Halos 

We consider a flat ACDM universe (Ha + f2 m =1), with the 
following values of the cosmological parameters: f2,„ = 0.27 
and H = 71 km s" 1 Mpc" 1 dSoergel et al.ll2003l) . In a flat uni- 
verse, the critical density can be written as 
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The virial radius of a DM halo with a virial mass mpM is then 
as follows: 
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where the spherical collapse overdensity, A c , and the matter 
density of the universe in units of critical density at the red- 
shift of z, $1%, are given by 



A e = 18tt 2 + 82^-39x 2 
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Here r s and c = r v \ r /r s are the scale radius and concentration 
of the halo. Burkert halos, on the other hand, have a flat core. 
Their outer density profile slope is identical to NFW halos 
(7 = -3). This model fits well the rotational cu rves of disk 
galaxies (Burkert 1995; Salucci & Burkert 2000) and has the 
following density profile: 



Pb(r) = 
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where 



Pom : 
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ln(l + c) + - ln(l + c ) - arctan c 
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The concentration c of cosmological halos has a weak de- 
pendence on a virial mass ( with less massive halos being 



more c oncentrated on average). Sternberg, McKee, & Wolfire 
(2002) give the following expression for a concentration of 
low mass DM halos in ACDM cosmological simulations at 
z = 0, which was obtained from the analysis of halos with the 
virial masses of 10 8 . . . 10 11 M Q : c = 27(m n ivi /10 9 M ( T 1 r a08 . 
Combining this expression with the result of B ullock et all 
(2001) that the concentration scales with a redshift as (l+z) -1 , 
we obtained the following formula for a concentration of low 
mass halos at different redshifts: 



OT DM 



27 / 

T+z ^ 10 9 M Q 
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(9) 



and 



iZhao et alJ l l2003albl) showed that CDM halos do not have 
concentrations smaller than ~ 3.5, hence equation be- 
comes invalid for halos with c < 4. Sternberg et al. (2002) 
demonstrated that equation describe s very well co ncentra- 
tions of four dwarf disk galaxies from Burkertl i ll 9951) fitted by 
a Burkert profile, so it can be used for both NFW and Burkert 
halos. 

2.3. Stellar Cores 

The initial non-equilibrium configuration of proto-GCs was 
modeled after MS04 as a homogeneous isothermal stellar 
sphere with isotropic Maxwellian distribution of stellar ve- 
locities. Stellar clusters of different mass initially had the 
same universal values of stellar density and velocity disper- 
sion: p it * = 14 Mq pc" 3 and cr fj * = 1.91 km s" 1 (MS04). As in 
MS04, we use a mass parameter (3 to describe different proto- 
GC models: 



l-0„ 



a„(i+z) 3 



(4) 



llBarkana & Loebl2001l) . Here x=fl z m -l. 

We consider two types of DM density profiles: NFW ha- 
los and Burkert halos. The NFW model describes reasonably 
well DM halos from cosmological CDM simulations. It has 
a cuspy inner density profile with a slope of 7 = -1, and a 
steeper than isothermal outer density profile: 



10V 



/ 375 



\47T/9; i *G 3 
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(10) 



where G is the gravitational constant. The connection be- 
tween the initial virial parameter v and (3 is v = 2K*/W* = 
10-W 3 , where and W* are initial kinetic and potential en- 
ergies of the stellar cluster. 

We use models from MS04 to plot in Figure^a dependence 
of the fraction of stars, which form a gravitationally bound 
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FIG. 1 . — Fraction of stars remaining gravitationally bound after the initial 
relaxation phase in proto-GC models of MS04. Empty circles mark models 
with /3 = -0.6 ("H"), p = 0.4 ("W"), and /3 = 1.4 ("C"). 



cluster after the initial relaxation phase of a proto-GC, on the 
mass parameter (3 (or the virial parameter v). Please note that 
this is the case of no DM (bare stellar cores). As one can 
see, there are three different regimes of the initial relaxation 
phase: 1) "Hot" collapse (-0.7 < /3 < -0.35 or 3 > v > 1.7) 
results in a substantial loss of stars which initial velocity was 
above the escape speed for the system; only the slowest mov- 
ing stars collapse to form a bound cluster. 2) "Warm" collapse 
(-0.35 < < 0.85 or 1.7 > v > 0.27) is very mild and non- 
violent, and results in virtually no escapers. 3) "Cold" col- 
lapse (J3 > 0.85 or v < 0.27) is violent, and leads to increas- 
ingly larger fraction of escapers for colder systems. In the last 
case, the nature of unbinding of stars is very different from the 
"hot" case. During a cold collapse, the radial gravitational po- 
tential is wildly fluctuating. At the intermediate stages of the 
relaxation, the central potential of the cluster becomes much 
deeper than the final, relaxed value. As a result, a fraction of 
stars are accelerated to speeds which will exceed the escape 
speed of the relaxed cluster. Analysis shows that most of es- 
capers are the stars that initially were predominantly in the 
outskirts of the cluster and moving in the outward direction. 

These three regimes of the initial relaxation of proto-GCs 
are of very general nature, and are not restricted to the GC for- 
mation scenario of MS04. To illustrate that, let us write down 
an expression for the virial ratio of a homogeneous isothermal 
sphere, which is applicable to both pre-starburst gas phase of 
a forming GC, and the initial stellar configuration after the 
star burst takes place: 



Here M and R are the mass and the radius of the system, and 
(V 2 ) is the mean-square speed of either gas molecules or stars. 
Assuming that newly-born stars have the same velocity dis- 
persion as star-forming gas, equation il Q suggests three fol- 
lowing physical scenarios resulting in our "hot", "warm", and 
"cold" stellar configurations (in all scenarios we start with an 
adiabatically contracting gas cloud of a Jeans mass, which 
corresponds to i/ gas = 1): 1) non-efficient stars formation with 



FIG. 2. — Redshift dependence of the minimum mass of a stellar core, 
with the universal initial density p; * = 14 Mq pc~ 3 , which dominates DM 
at the center of a DM halo. Thick lines (colored green in electronic edition) 
correspond to stellar cores in NFW halos, whereas thin lines (colored blue 
in electronic edition) correspond to clusters in Burkert halos. Long-dashed, 
solid, and short-dashed lines correspond to \ = 0.0025,0.0088,0.2. 



a subsequent loss of the remaining gas leads to our "hot" case, 
2) almost 100%-efficient star formation results in the "warm" 
case, and 3) runaway cooling of the whole cloud on the time 
scale shorter than the free-fall time, leading to high-efficiency 
star formation, results in our "cold" case. 

2.4. Physical Parameters of the Models 

In this paper, we consider three following initial stellar 
configurations, which can be considered to be typical "hot", 
"warm", and "cold" cases (see Figure 0: (3 = -0.6 (model 
"H", stands for "Hot"), j3 = 0.4 (model "W", stands for 
"Warm"), and = 1.4 (model "C", stands for "Cold"). We 
study relaxation of stellar cores in either NFW (suffix "«"; 
e.g., "W„") or Burkert (suffix "b") DM halos, which makes 
up a total of 6 different combinations of stellar cores and DM 
halos. The physical parameters of these 6 models are summa- 
rized in Table [I] 

The stellar masses in our models cover the whole range 
of GC masses (see Table The initial stellar radius r* is 
much smaller than the scale radius of the DM halo r s in all 
the models. The masses of DM h alos mpM range fro m 10 6 
to 10 s M Q . From Figure 6 of Barkana & Loeb (2001), these 
DM halos correspond to 1 - 1 .5(7 fluctuations collapsing at the 
redshift of z = 7. As one can see, halos in this mass range are 
populous enough at z ~ 7 to account for all observed GCs. 

In all our models, stars initially dominate DM in the core 
(S = ffit/mDM^] > 1, see Table[U, which can be considered 
as a desirable property for a star formation to take place. It is 
instructive to check if the above condition (S > 1) holds for 
other plausible values of the GC formation redshift z and stel- 
lar mass fraction \- I n Figure|2]we show the redshift depen- 
dence of the minimum mass of a stellar core in a DM halo, 
satisfying the condition S > 1, for a few fixed values of \ 
— 0.0025, 0.0088 (our fiducial value), and 0.20. Both NFW 
(thick lines) and Burkert (thin lines) cases are shown. As one 
can see, only in the extreme case of an NFW halo at z = 20 
with x = 0.0025 fraction of total mass in stars, the minimum 
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TABLE 1 

Physical parameters of the models 



Model 


/3 




m. 


r* 


T* 


' *.min 


»1DM 


c 


r vir ?s 


Of.DM 


TJM 


S 








M 


pc 


Myr 


pc 






pc pc 


pc 


Myr 




W„ 


0.4 


0.54 8 


1.8 x 10 4 


11.2 


0.49 


2.98 


io 7 


4.88 


885 181 


395 


52 


4.7 


W b 


0.4 


0.54 8 


;.8 x io 4 


11.2 


0.49 


2.98 


10 7 


4.88 


885 181 


436 


61 


117 


c„ 


1.4 


0.12 8 


;.8 x io 5 


24.2 


0.078 


1.57 


IO 8 


4.06 


1906 470 


892 


56 


5.8 


Q, 


1.4 


0.12 8 


;.8 x io 5 


24.2 


0.078 


1.57 


10 s 


4.06 


1906 470 


997 


66 


173 


H„ 


-0.6 


2.5 8 


;.8 x 10 3 5.21 


17 


4.14 


IO 6 


5.87 


411 70 


174 


48 


3.7 


H fc 


-0.6 


2.5 8 


i.8 x 10 3 


5.21 


17 


4.14 


10 s 


5.87 


411 70 


190 


55 


78 



NOTE. — Here v* and r* are the initial virial ratio (assuming, that there is no DM) and radius of the stellar core; r* is the crossing time at the half-mass 
radius for the relaxed, purely stellar models from MS04 rescaled to p,> = 14 Mq pc~ 3 and cr, > = 1.91 km s _1 ; r* ra j n is the minimum attained half-mass radius for 
relaxing stellar clusters from MS04 (rescaled to the above values of p, „ and cr, r/, dm and "TDM are the half-mass radius and the crossing time at the half-mass 
radius for DM halos; S = m* /mdm^*) is the ratio of the stellar mass to the mass of DM within the initial stellar radius r„. 



mass of a dominant stellar core m* m in = 1.7 x 10 3 M Q is ap- 
proaching the mass range of GCs. The conclusion we arrive 
at is that, for all plausible values of z and x< stellar cores with 
the universal initial density p,- * = 14 Mq pc" 3 and GC masses 
dominate DM at the center of DM halos. 

2.5. Setting Up N-Body Models 

DM halos in our models are assumed to have an isotropic 
velocity dispersion tensor. To generate A^-body realizations of 
our DM halos, we sample the following probability density 
function (PDF) JWidrowl2000l) : 

P(E,R)<xR 2 (^-E) l/2 F(E). (12) 

Here ^ and E are the relative dimensionless potential and 
particle energy in units of AnGpa. n r 2 (for NFW halos) and 
TT 2 Gpoj,r 2 (for Burkert halos), F(E) is the phase-space dis- 
tribution function, and R = r/r s is the dimensionless radial 
coordinate of the particle. The dimensionless potential ^ 
is equal 1 at the center and in infinity. For NFW halos, 
'I' = R' 1 ln(l + R). For Burkert halos, the corresponding ex- 
pression is given in Appendix lAl (eg. IA1| ). We use the an- 
alytic fitting formulae of iWidrowl (120001) to calculate F(E) 
for NFW model, and derive our own fitting formulae for the 
Burkert profile (see AppendixlAV Our models are truncated at 
the virial r adiu s r wu (eq. 0): F(E) = for r > r v i r . As we will 
show in § 13.21 this truncation results in an evolution in the 
outer DM density profile (with the profile becoming steeper 
for r ^> r s ), which should not affect our results. 

We sample the PDF (eq. 11121 ) in two steps: 1) A uniform 
random number x € [0 ... 1] is generated. The radial distance 
R is then obtained by solving numerically one of the two fol- 
lowing non-linear equations: 

\n(l+R)-R/(l+R) = x[ln(l+c)-c/(l+c)] (13) 
for the NFW profile and 

ln( 1 + R) + [ln( l+R 2 )]/2- arctanfl = 

= x [ln( 1 + c) + [ln( 1 + c 2 )] /2 - arctanc] (14) 

for the Burkert profile. 2). Given R, we can calculate 'J7. 
Instead of equation (I12> . we can now sample the following 
PDF: 

fty-E\ 1/2 F(E) 
P(E)(xU={ — — 1. (15) 



As F(E) is a monotonically increasing function and < E < 
'J?, the following inequality holds: < II < 1. Two uniform 
random numbers, y 6 [0 . . . 1] and E € [0 . . . $], are generated. 
If y < IT, we accept the value of the dimensionless energy E; 
otherwise, we go back to the previous step (generation of y 
and E). 

The velocity module v can be obtained from the equation 
for the total energy of a particle: E = v 2 /2 + ^. Finally, two 
spherical coordinates angles 8 and ip are generated for both 
radius and velocity vectors using the following expressions: 
cos 8 = 2t-l and ip = 2iru. (Here t and u are random numbers 
distributed uniformly between and 1.) 

The above method to generate TY-body realizations of ei- 
ther NFW or Burkert halos does not use a "local Maxwellian 
approximation" to assign velocities to particles. Instead, it 
explicitly uses phase-space distribution functions, which was 
shown to be a superior way of se tting up A^-body models 
( Kazantzidis. Magorrian. & Moore 2004). 

Stellar cores are set up as homogeneous spheres at the cen- 
ter of DM halos: R = r*x 1//3 , where x £ [0. . . 1] is a uniformly 
distributed random number. We use equal mass stellar parti- 
cles. Velocities of the particles have a Maxwellian distribu- 
tion. The components of the radius and velocity vectors are 
generated in the same fashion as for DM particles. 

2.6. Numerical Parameters of the Models 

In addition to 6 models from Table ^ which consist of a 
live DM halo plus a stellar core, we also ran 2 more models 
with a static DM potential (marked with the letter "S" at the 
end). The numerical parameters for these 8 models are listed 
in Table [2] In this table we also list parameters for the three 
properly rescaled stars-only models from MS04 (Ho, Wo, and 
Co), which will be used as reference cases. 

To run our models, we use a pa rallel version of the multi- 
stepping tree code GADGET- 1.1 ( Springel et alJl2001l) . The 
code allows us to handle separately stellar and DM particles, 
with each species having a different softening length — e* and 
cdm, respectively. We slightly modified the code by introduc- 
ing an optional static DM potential (either NFW, or Burkert). 

The values of the softening lengths were chosen to be com- 
parable with the initial average interparticle distance. We 
used the following expression for the stellar particles: e = 
0.77 r h N-^ 3 iHavashiet ail 12003). where the half- mass ra- 
dius r;, was measured at the initial moment of time (so r/, = 
r^/2 1 ^). The same expression was used to calculate €dm 

for 

H„ t b models, whereas for the models W„ and C„j, we used 
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TABLE 2 

Numerical parameters of the models 



Model 


N* 


e* 






t] 


h 


Af raax 


Note 






pc 




pc 


Myr 


Myr 


Myr 




W 


10 5 


0.15 






190 


370 


0.3 


1 


W„S 


10 4 


0.30 






95 


190 


0.2 


2 


w„ 


It) 4 


0.30 


10 6 


1.5 


200 


250 


0.2 






10 4 


0.30 






95 


190 


0.2 


2 


w b 


10 4 


0.30 


10 6 


1.5 


150 


250 


0.2 




Co 


10 5 


0.32 






56 


92 


0.3 


1 


c„ 


10 s 


0.32 


5 x 10 5 


4.7 


130 


190 


0.2 




c b 


10 5 


0.32 


5 x 10 5 


5.2 


95 


190 


0.2 




Ho 


10 5 


0.069 






1300 


1800 


0.3 


1 


H„ 


1() 4 


0.15 


5 x 10 5 


1.7 


1500 


3700 


20 




H fc 


10 4 


0.15 


5 x 10 5 


1.8 


1500 


3700 


20 





NOTE. — 1) Models from MS04 rescaled to p,.„ = 14 Mq pc 3 and 
o-,.* = 1.91 km s~' (stars only, no DM). 2) Static DM potential. Here A/» 
and /Vdm are the number of stellar and DM particles; e* and eou are the 
softening lengths for stars and DM; t\ , t2, and A? max are the moment of time 
when most of model parameters converge to their quasi-steady-state values, 
the total evolution time, and the maximum value for individual timesteps. 



half of the value given by the equation of lHavashi et al J (120031) 
to achieve a better spatial reso lution within the stellar core. 

Aarseth. Lin. & Papaloizou ( 1988) showed that collapsing 
homogeneous stellar spheres with a non-negligible velocity 
dispersion er,- * do not experience a fragmentation instability 
(and as a consequence preserve orbital angular momentum of 
stars), if the following condition is met: /V* > [m*/(r*of*)] 3 , 
where 7V» is the number of stars. For our models, this adia- 
baticity criterion can be rewritten as 

logiV* >2/3+31og5. (16) 

As we discussed in MS04, if real GCs indeed started off 
as isothermal homogeneous stellar spheres, they should have 
collapsed adiabatically. We also showed that if both the con- 
dition in equation dl6l and e» < 0.25r*, m i n are met, the results 
of simulations of collapsing stellar spheres do not depend on 
the number of particles Af* and the softening length e*. (Here 
r*,min is the minimum half-mass radius of the cluster during 
the collapse.) By comparing Tables ^ and |2] one can see that 
the values of N* and e* we use in our simulations meet both 
above conditions. 

The evolution time ?2 in our runs is at least 3 times longer 
than the crossing time for DM tdm, and hundreds and even 
thousands times longer than the crossing time for stellar par- 
ticles r» (see Tables [2 and |3- As we will see in the next sec- 
tion, this time is enough for stellar and DM density profiles 
to converge. On the other hand, it is short enough to avoid 
significant dynamic evolution in the stellar clusters caused by 
encounters between individual particles. 

The individual time steps in the simulations are equal to 
(2?ye/a) 1 ' /2 , where a is the acceleration of a particle, and pa- 
rameter r\ controls the accuracy of integration. We used a very 
conservative value of r\ = 0.0025 and set the maximum possi- 
ble individual time step Af max value to either 20 or 0.2 Myr 
(see Table|2}. As a result, the accuracy of integration was very 
high: AE tot < 0.06%, where AE tot is the maximum deviation 
of the total energy of the system from its initial value. One 
has to keep in mind though that numerical artifacts are the 
most pronounced in the central densest area, where the fre- 



quency of strong gravitational interaction between particles 
is the highest. We estimate the severity of these effects by 
looking at the total energy conservation in our purely stellar 
models from MS04: in models Ho, Wo, and Co the values of 
AEtot are - 0.05%, - 0.85%, and ~ 0.8% (for the model W 
we used a larger value of 77 = 0.02, hence the relatively large 
errors). As you can see, even though the numerical artifacts 
in the dense stellar cluster are more visible than in the DM + 
a stellar core case, the magnitude of these effects is still rea- 
sonably low. 

For completeness sake, we also give the values of other 
code parameters which control the accuracy of simula- 
tions: ErrTolTheta=0.6, TypeOfOpeningCriterion=l, Er- 
rTolForceAcc=0.01, MaxNodeMove=0.05, TreeUpdateFre- 
quency=0.1, and DomainUpdateFrequency=0.2. (Please see 
the code manual 1 for explanation of these parameters). The 
code was compiled with the option DBMAX enabled, which 
allowed it to use a very conservative node-opening criterion. 

3. RESULTS OF SIMULATIONS 
3.1. General Remarks 

For each of our models, we generated 100-200 individual 
timeframe snapshots. The last 20-60% of the snapshots (cor- 
responding to the time interval ti ... ti, see Table|2j were used 
to calculate the properties of relaxed models, including radial 
profiles and different stellar cluster parameters listed in Ta- 
ble[3] In all snapshots, we removed in an iterative manner the 
escapers (particles whose velocity exceeds the escape velocity 
Ksc = [ _ 2$] l / 2 , where $ is the local value of the gravitational 
potential). This procedure affected models from MS04 (Ho, 
Wo, and Co), and DM particles only in the models presented 
here. One of the reasons for removing escapers was to make 
the results of the simulations directly comparable with the re- 
sults from Paper II where this procedure is applied to discount 
particles stripped by tidal forces. 

The radial profiles shown in this section were obtained by 
averaging the corresponding profiles from individual snap- 
shots in the time interval t\ . . .to. For each profile, we only 
show the part which is sufficiently converged (the dispersion 
between different snapshots is very small). For projected 
quantities (such as surface brightness £y and line-of-sight ve- 
locity dispersion o>) we use the projection method described 
in Appendix 151 which produces much less noisy radial pro- 
files than a straightforward projection onto one plane or three 
orthogonal planes. This is crucial for the analysis of the prop- 
erties of the outer, low density parts of stellar clusters, which 
is the main purpose of this study. 

The virial ratio v for our models is consistent with unity 
within the measurement errors at the end of simulations. All 
the global model parameters (listed in Table [5} converge to 
their final values by t = ?2- Similarly, radial density and veloc- 
ity dispersion profiles converge to their final form over a wide 
range of radial distances. We explicitly checked that the evo- 
lution time t2 is more than three crossing times at the largest 
radius for stellar density profiles shown in Figures|3]|5] and0 
All the above let us conclude that the results presented here 
are collisionless steady-state solutions. 

3.2. Warm Collapse 

In the absence of DM, an isothermal homogeneous stellar 
sphere with the mass parameter (3 = 0.4 relaxes to its equi- 
librium state with virtually no stars lost (see Figure [Q. This 

' |http : // www . mpa-gar ching .mpg. de/ gadget /| 
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TABLE 3 

Parameters for relaxed stellar clusters 



Model 


Th,* 




Pc 


R/i,* Rhb 


O"0 


XV 


'"0 


/o 


T 


r P 


r„, 




pc 


km s _1 


M Q pc 3 


pc pc 


km s 


mag arcsec~ 2 


pc 




M L„' 

W (.J 


pc 


pc 


W 


4.64 


4.07 


310 


3.64 2.77 


3.84 


18.71 


3.01 


0.226 


1.43 






W„S 


4.51 


4.36 


360 


3.50 2.63 


4.11 


18.63 


2.96 


0.253 


1.60 






w„ 


4.08 


4.83 


500 


3.16 2.31 


4.55 


18.39 


2.78 


0.279 


1.80 


7.74 


19.7 




4.67 


4.05 


310 


3.64 2.81 


3.83 


18.73 


2.98 


0.238 


1.43 






Wj, 


4.65 


4.12 


350 


3.61 2.54 


3.89 


18.65 


2.86 


0.229 


1.52 


17.8 


54.7 


Co 


2.55 


14.4 


1.3 X 10 4 


2.03 1.52 


13.6 


15.31 


1.64 


0.241 


1.43 






C„ 


4.50 


14.3 


1.1 x 10 4 


3.54 1.55 


13.4 


15.45 


1.75 


0.166 


1.55 


10.0 44.0 


c h 


4.50 


14.4 


1.3 x 10 4 3.52 1.49 


13.5 


15.30 


1.62 


0.168 


1.42 


27.6 


125 


Ho 


20.2 


0.54 


0.71 


14.9 7.13 


0.51 


24.28 


8.29 


0.154 


1.64 






H„ 


5.47 


1.57 


25 


4.20 2.79 


1.56 


21.44 


4.06 


0.329 


2.88 






Hi, 


11.5 


0.97 


4.8 


8.84 4.43 


0.97 


22.70 


5.71 


0.212 


2.23 


10.4 


18.2 



NOTE. — Here r;, », a c , and p c are half-mass radius, central velocity dispersion, and central density of stellar clusters; /?/,,„, Rhb, and XV are projected 
half-mass radius, half-brightness radius (where surface brightness is equal to one half of the central surface brightness), projected central dispersion, and central 
surface brightness (see eq. 1171 ): rp and fo are the King core radius ro = [9cr 2 /(4-irGpc)] 1 / 1 and the fraction of the total stellar mass inside ro; T is the apparent 
central mass-to-light ratio (see eq. 1181 ); r p is the radius where density of DM and stars becomes equal; r,„ is the radius where enclosed mass of DM becomes 
equal to that of stars. 




1 2 3 1 2 

log (r), pc log (r), pc 



FIG. 3. — Radial density profiles for warm collapse models, (a) The case of NFW halo, (b) The case of Burkert halo. Thick lines correspond to DM density; 
thin lines (colored red in electronic edition) show stellar density. Solid lines correspond to the cases of a live DM halo + a stellar core (models W„ and Wft), 
long-dashed lines depict analytic DM profiles and a relaxed stellar cluster profile in the absence of DM (model Wo), and dotted lines show stellar density profiles 
for static DM models (W„S and W^S). Vertical long-dashed, short-dashed, and dotted lines mark the values of edm, f m , and r s . 



makes it an interesting case to test the result o flPeeblesI i ll 9841) . 
that a stellar cluster inside a static constant density DM halo 
can acquire a radial density cutoff similar to that expected 
from the action of tidal forces of the host galaxy, for the case 
of live DM halos with cosmologically relevant density profiles 
and initially non-equilibrium stellar cores. 

In Figure 13 we show the radial density profiles for relaxed 
stellar core and DM halo, both for NFW model (panel a) and 
Burkert model (panel b), for the case of a warm collapse ("W" 
models). Here long-dashed lines show the density profiles for 
stars in the absence of DM (thin lines) and DM in the absence 
of stars (thick lines). Solid lines, on the other hand, show 
the density profiles for stars and DM having been evolved to- 
gether. We also show as dotted lines the stellar radial density 
profiles for the case of a static DM potential. 



As you can see from Figure|5] in the absence of DM (model 
Wo), the relaxed stellar cluster has a central flat core and a 
close to power-law outer density profile. The radial density 
profile has also a small "dent" outside of the core. It is not 
a transient feature (in the collisionless approximation), as the 
virial ratio for the whole cluster is equal to unity within the 
measurements errors (y = 1 .001 ±0.005 for t = t\ . . J2) and the 
evolution time is > 7 crossing times even at the very last radial 
density profile point. As we show in Paper II, collisional long- 
term dynamic evolution of a cluster will remove this feature, 
bringing the overall appearance of the density profile closer to 
that of the majority of observed GCs. The central stellar den- 
sity in the model Wo is ~ three orders of magnitude larger 
than the central density of DM in the undisturbed Burkert 
halo, and becomes comparable to the density of the undis- 
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FIG. 4. — Radial profiles of observable quantities for warm collapse models, (a) Surface brightness Ey. (b) Line-of-sight velocity dispersion o>. Thick 
lines (colored green in electronic edition) correspond to NFW cases; thin lines (colored blue in electronic edition) correspond to Burkert cases. Solid lines show 
profiles for models W„ and short-dashed lines correspond to models W„S and WjS, and long-dashed lines correspond to model Wo. Vertical short-dashed 
lines mark the values of r„, . 



turbed NFW halo only at very small radii r < 0. 1 pc (which is 
outside of the range of radial distances in Figure^). 

In the case of a live DM halo coevolving with the stellar 
core, the DM density profile is significantly modified within 
the central area dominated by stars (see Figure [3}. In both 
NFW and Burkert halos, DM is adiabatically compressed by 
the potential of the collapsing stellar cluster to form a steeper 
slope in the DM radial density profile. The slope of the inner- 
most part of the DM density profile is 7 ~ —1 .5 for the NFW 
halo and 7 ~ -1 .0 for the Burkert halo. The lack of resolution 
does not allow us to check if the slope stays the same closer 
to the center of the DM halos. 

In the outer parts of the DM halos (beyond the scale radius 
r s , which is marked by vertical dotted lines in Figure [3), the 
DM density profile becomes steeper than the original slope of 
7 = -3. This can be explained by the fact that our DM halo 
models are truncated at a finite radius r wu and hence are not in 
equilibrium. It should bear no effect on our results, as in all 
our models stellar clusters do not extend beyond r s . 

As you can see in Figure[5] the impact of the presence of a 
live DM halo on the stellar density profile is remarkably small 
(especially for the Burkert halo). In the presence of DM, the 
central stellar density becomes somewhat larger (by ~ 60% 
and ~ 10% for the NFW and Burkert halos, respectively; see 
Tabled. 

The most interesting effect is observed in the outer parts of 
the stellar density profiles, where the slope of the profiles be- 
comes significantly larger than th at of a purely st ellar cluster, 
which is similar to the result of Peebles (1984). The den- 
sity profile starts deviating from the profile for the model Wo 
somewhere between the radius r p and the radius r m (where the 
enclosed masses for stars and DM become equal, see Table[3}. 

There is a simple explanation for the steepening of the stel- 
lar density profile seen in Figure [3] around the radius r m . In 
the absence of DM, warm collapse of a homogeneous stellar 
sphere (our model Wo) is violent enough to eject a number of 
stars into almost radial orbits, forming a halo of the relaxed 
cluster. In the case when the DM halo is present, dynamics 



does not change much near the center of the cluster where 
stars are the dominant mass component. Around the radius r m 
the enclosed DM mass becomes comparable to that of stars. 
As a result, gravitational deceleration, experienced by stars 
ejected beyond r m , more than doubles, resulting in a signif- 
icantly smaller number of stars populating the outer stellar 
halo at r > r m . This should lead to significantly steeper radial 
density profile beyond r m , as observed in our W„ j, models. 

As in lPeeblesl tl984) model, in our model the steepening of 
the stellar density profile is caused by the fact that at large 
radii the gravitational potential is dominated by DM. The 
princip al difference between the two models is that in Peebles 
( 1984) the stellar cluster is assumed to be isothermal and in 
equilibrium; our clusters are not in equilibrium initially, and 
their final equilibrium configuration is not isothermal (stars 
are dynamically colder in the outskirts of the cluster). 

To allow a direct comparison of the model results with ob- 
served GCs, in Figure 0] we plot both the surface brightness 
S v profile (panel a) and the line-of-sight velocity dispersion 
<r r profile (panel b) for warm collapse models. The surface 
brightness Ey in units of mag arcsec" 2 is calculated using the 
following formula: 

E l , = y -2.5[log(C/T GC ) + 2-21og(36OOx 180/tt)]. (17) 

Here V Q = 4™87 is the absolute V-band magnitude of the 
Sun, £ is the projected surface mass density in units of 
Mq pc" 2 , and Tec is the assumed V-band mass-to-light ratio 
for baryons in GCs. For Tqc we a dopt the observation allv de- 
rived averaged value of 1.45 from McLauahlin ( 1999), which 
is also between the two stellar-synthesis model values (1.36 
for the Salpeter and 1.56 for the composite — with a zero 
sl ope for stellar mass es < 0.3 Mq — initial mass functions) 
of lMateo et aT](f7998b . 

As you can see in Figure |4] the presence of a DM halo in- 
troduces an outer cutoff in the surface brightness radial profile 
of stellar cluster, which makes the profile look very similar to 
a King model profile. Interestingly, stellar cluster evolving in 
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a live NFW DM halo acquires even steeper outer density and 
brightness profiles than for the case of a static DM potential 
(see Figures|5t and|4^). 

Analysis of the line-of-sight velocity dispersion profiles for 
warm collapse models (Figure^) shows that the presence of 
a DM halo slightly inflates the value of a> in the core of the 
cluster (by < 20%, see Table[3j. The outer dispersion profile 
either stays unchanged (as in the case of the Burkert halo), 
or becomes slightly steeper (for the NFW model). The static 
DM models present profiles of an intermediate type. 

A slight increase of o> in the central part of the cluster can 
be misinterpreted observationally as a GC with no DM which 
has a bit larger value of the baryonic mass-to-light ratio T. To 
quantify this effect, we apply a core fitting (or King's) method 
of finding the central value o f T in spherical stellar systems 
dRichstone & Tremainell9 86): 



InGhRhb 

Here To is the central surface brightness in physical units. We 
assume here that Io = (o/Tgc, where Q) is the projected sur- 
face mass density at the center of the model cluster. The core 
fitting method assumes that the mass-to-light ratio is inde- 
pendent of radius, which is obviously not true for our stars + 
DM models. In our models the stellar particles have the same 
mass, so there is no radial mass segregation between high and 
low mass stars caused by the long-term dynamic evolution of 
the cluster. As a result, the equation (II 81 should not be used 
to compare our models with real GCs. Instead, we use it to 
check if there is a change in the apparent mass-to-light ratio T 
between our purely stellar models and the models containing 
DM. We list the values of T for different models in Table [3] 
As you can see, in the case of NFW halo (model W w ), the 
presence of the DM halo leads to ~ 25% increase in the value 
of the apparent mass-to-light ratio. Such a small increase is 
well withi n the observed d ispersion of T values for Galactic 
GCs jPrvor & Mevlarl 19931) . In the case of the Burkert halo, 
the apparent central mass-to-light ratio is only ~ 6% larger 
than for the purely stellar model. 

3.3. Cold Collapse 

A principal difference between warm and cold collapse 
models from MS04 is that in the cold collapse case a signifi- 
cant fraction of stars becomes unbound after the initial relax- 
ation phase (see Figure [Q. In our model Co with the mass 
parameter (3 = 1 .4 the fraction of the stars lost is 30%. The ra- 
dial density profile for Co model is similar to that of Wo model 
(see Figure |5}- One can intuitively expect the escapers from 
the model Co to be trapped inside our models containing DM, 
forming a distinctive feature in the outer density and velocity 
dispersion profiles. 

As you can see in Figures[3]and[6] our cold collapse models 
do show such features. Both density and surface brightness 
profiles become more shallow in the outer parts of the stellar 
clusters where DM dominates stars. This is valid for both 
NFW and Burkert halos. For the NFW and Burkert cases, 
the slope of the outer stellar density profile is 7 = -2.6 and 
-2.2, respectively. For the purely stellar case (model Co) the 
slope is 7 = -3.8, so the relative change in the slope is A7 = 
1.2 and 1.6 for the NFW and Burkert cases. This behavior 
is mimicked by the corresponding surface brightness profiles 
(Figure|6t). It is interesting to note that, for the C„ model the 
radial Y,y profile exhibits a significant steepening of the slope 



in the outmost parts of the cluster, creating an appearance of 
a "tidal cutoff" (similarly to the warm collapse cases). 

Even more pronounced are features in the line-of-sight ve- 
locity dispersion profiles (Figure^). In both NFW and Burk- 
ert cases, there is a plateau in the radial <r, profiles around 
the radius r m where DM becomes the dominant mass compo- 
nent. The apparent "break" in the surface brightness profile 
and accompanying flattening of the radial a r profile, seen in 
our models containing DM, can be misinterpreted observa- 
tionally as a presence of "extratidal" stars heated by the tidal 
field of the host galaxy. 

Similarly to warm collapse models, the apparent mass-to- 
light ratio T for cold collapse models in the presence of a live 
DM halo is very close to the purely stellar case Co (see Ta- 
ble [3}- Interestingly, the half-mass radii for warm and cold 
collapse models with DM are almost identical. This is con- 
sistent with the properties of the observed GCs, which have 
comparable half-mass radii for a wide range of cluster masses. 

DM density profiles for the cold collapse cases (Figure [5} 
exhibit a very similar behavior to the warm collapse models in 
the stars dominated central region. For both NFW and Burk- 
ert halos, the innermost slopes of the density profiles become 
steeper (7 ~ -1.8 and -0.6, respectively) in the presence of 
a stellar core. The NFW DM density profile shows a break 
around the radius r m . 

As you can see in Tableland Figure|5J the final stellar clus- 
ter half-mass radius is smaller than the DM softening length 
£dm ~ 5 pc. We reran models C„ and Q, with much smaller 
value of 6dm = 1 pc to test the possibility that our results were 
influenced by the fact that DM is not resolved on the stellar 
cluster scale. For both C„ and Q,, the relaxed stellar profiles 
are found to be virtually identical to the cases with larger 6dm- 
In particular, a "kink" seen in Figures [5^ and 13) around the 
radius cdm is also present at the same location in the simula- 
tions with much smaller value of cdm, and is definitely not a 
numerical artifact. In the case of NFW halo, the central stel- 
lar velocity dispersion becomes slightly larger, which results 
in somewhat larger value of T = 1.88. This could be in part 
because of artificial mass segregation, which should be more 
pronounced in the case of 6dm being significantly lower than 
the optimal value (in our cold collapse models, DM particles 
are ~ 20 times more massive than stellar particles). As a re- 
sult, the actual T value could be even closer to the baryonic 
value. In the case of Burkert halo with £dm = 1 pc, the core 
mass-to-light ratio T = 1.45, which is identical to the case of 
no DM. We conclude that our choice of £dm did not affect the 
main results presented in this Section. 

3.4. Hot Collapse 

It is generally assumed that in an equilibrium star-forming 
cloud no bound stellar cluster will be formed after stellar 
winds and supernova explosions expel the remaining gas if 
the star formation efficiency is less than 50% (which corre- 
sponds to the initial virial parameter for the stellar cluster of 
v > 2). In MS04 we showed that it is not true for initially 
homogeneous stellar clusters which have a Maxwellian dis- 
tribution of stellar velocities. In such clusters with the initial 
virial parameter as large as 2.9 (corresponding to the mass 
parameter (3 as low as -0.7), a bound cluster, containing less 
than 100% of the total mass, is formed by th e slowest moving 
stars. (A similar conclusion was reached by Boilv & Kroupa 
2003 for initially polytropic stellar spheres.) 

In our model Ho (J3 = -0.6, v = 2.5), ~ 57% of all stars stay 
gravitationally bound after the initial relaxation phase (see 
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FIG. 5. — Radial density profiles for cold collapse models, (a) The case of NFW halo, (b) The case of Burkert halo. Thick lines correspond to DM density; 
thin lines (colored red in electronic edition) show stellar density. Solid lines correspond to the cases of a live DM halo + a stellar core (models C„ and C/,); 
long-dashed lines depict analytic DM profiles and a relaxed stellar cluster profile in the absence of DM (model Co). Vertical long-dashed, short-dashed, and 
dotted lines mark the values of edm, I'm, and r s . 




FIG. 6. — Radial profiles of observable quantities for cold collapse models, (a) Surface brightness Xy. (b) Line-of-sight velocity dispersion o>. Thick lines 
(colored green in electronic edition) correspond to NFW cases; thin lines (colored blue in electronic edition) correspond to Burkert cases. Solid lines show 
profiles for models C„ and Q,; long-dashed lines correspond to model Co. Vertical short-dashed lines mark the values of r m . 



Figure 0. Because in both cold and hot collapse models for 
purely stellar clusters a significant fraction of stars become 
unbound, one might expect that stellar density and velocity 
dispersion profiles to be similar for these two cases in the pres- 
ence of DM. Analysis of Figures0and|S]shows that it is defi- 
nitely not the case: H n ^ models show completely different be- 
havior from C„,i models. The explanation is simple: Co model 
starts with the stellar density which is already larger than the 
density of DM within the stellar core in models C,,,/, (in other 
words, S > 1, see Table and after the relaxation phase 
reaches even higher density which is ~ (r*/r/, ») 3 /2 ~ 430 
times higher than the initial density. The hot collapse model 
Hq, on the other hand, having started with S > 1, acquires a 



density which is ~ 120 times lower than the initial density 
after the expansion and relaxation phase, resulting in S < 1 
for the relaxed stellar cluster even for the Burkert halo. Un- 
der these circumstances, radial profiles for H„ and H& models 
should be significantly different from the case with no DM 
(Ho), which is what is seen in FiguresQand[8] The only com- 
mon feature between the models C„,z, and H„ j, is that in both 
cases all would-be escapers are trapped by the potential of the 
DM halo. 

In the case of the NFW halo, central stellar density and cen- 
tral velocity dispersion become larger in the presence of DM 
by factors of 35 and 3, respectively. In the Burkert halo case, 
the increase is not as dramatic, but still significant: central 
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FIG. 7. — Radial density profiles for hot collapse models, (a) The case of NFW halo, (b) The case of Burkert halo. Thick lines correspond to DM density; 
thin lines (colored red in electronic edition) show stellar density. Solid lines correspond to the cases of a live DM halo + a stellar core (models H„ and H/,); 
long-dashed lines depict analytic DM profiles and a relaxed stellar cluster profile in the absence of DM (model Ho). Vertical long-dashed, short-dashed, and 
dotted lines mark the values of cdm, r„, (only the Burkert case is shown, as in the NFW case the enclosed DM mass is larger than that of stars at any radius), and 




FIG. 8. — Radial profiles of observable quantities for hot collapse models, (a) Surface brightness Ey, (b) Line-of-sight velocity dispersion <r, . Thick lines 
(colored green in electronic edition) correspond to NFW cases; thin lines (colored blue in electronic edition) correspond to Burkert cases. Solid lines show 
profiles for models H„ and long-dashed lines correspond to model Ho. Vertical short-dashed line marks the value of r„, (only the Burkert case is shown, as in 
the NFW case the enclosed DM mass is larger than that of stars at any radius); vertical dotted line corresponds to r s . 



density and velocity dispersion become 7 and 2 times larger 
than those in the absence of DM. Interestingly, as in the cold 
collapse cases, the presence of DM brings the half-mass ra- 
dius of the stellar core in H„ ^ models closer to the half-mass 
radius of the warm collapse models (see Table[3}. The appar- 
ent central mass-to-light ratio T is noticeably inflated by the 
presence of DM (especially for the case of NFW halo) — by 
76% and 36% for NFW and Burkert models, respectively. 

As you can see in Figure^, in the model H„ stars have a 
comparable density to DM in the core, and become less dense 
than DM outside of the core radius r^. The situation is dif- 
ferent (and more in line with warm and cold collapse mod- 



els) for the Hi, model (Figure |7J>), where stars are the dom- 
inant (though not by a large margin) mass component in the 
core. There are no obvious features in the density and surface 
brightness profiles caused by the presence of DM: the profiles 
look very similar to Wo and Co models which do not have 
DM. 

Similarly to warm and cold collapse cases, the presence of a 
stellar cluster makes DM denser in the central area. In the case 
of NFW halo, the innermost slope of the DM density profile 
is comparable to the slope of 7 = -1 for undisturbed halo. 
For the Burkert halo, the slope becomes steeper, reaching the 
value of 7 ~ -0.4 at the innermost resolved point. 
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The most interesting behavior is exhibited by the line-of- 
sight velocity dispersion profiles (Figure|SJ)). In the presence 
of DM, the profiles become remarkably flat, changing by mere 
0.2 - 0.3 dex over the whole range of radial distances (out to 
~ 15 apparent half-mass radii). The combination of the low 
mass, low central surface brightness, large apparent mass-to- 
light ratio, and almost flat velocity dispersion profiles seen in 
models H„ /, can be mistakingly interpreted as a GC observed 
at the final stage of its disruption by the tidal forces of the host 
galaxy. 

4. DISCUSSION AND CONCLUSIONS 

A significant hindrance to a wider acceptance of the primor- 
dial scenarios for GC origin is an apparent absence of DM in 
Galactic GCs. Many observational facts have been suggested 
to be evidence for GCs having no DM, including the presence 
of such features in the outer parts of the GC density profiles 
as apparent tidal cutoffs or breaks, relatively low values of 
the apparent central mass-to-light ratio T which are consis- 
tent with purely baryonic clusters, flat radial distribution of 
the line-of-sight velocity dispersion in the outskirts of GCs 
believed to be a sign of tidal heating, and non-spherical shape 
of the clusters in their outer parts. 

Here we present the results of simulations of stellar clus- 
ters relaxing inside live DM minihalos in the early universe 
(at z = 7). We study three distinctly different cases which can 
correspond to very different gas-dynamic processes forming 
a GC: a mild warm collapse, a violent cold collapse resulting 
in a much denser cluster with a significant fraction of stars 
escaping the GC in the absence of DM, and a hot collapse 
resulting in a lower density cluster with many would-be stars- 
escapers. We show that GCs forming in DM minihalos exhibit 
the same properties as one would expect from the action of the 
tidal field of the host galaxy on a purely stellar cluster: King- 
like radial density cutoffs (for the case of a warm collapse), 
and breaks in the outer parts of the density profile accompa- 
nied by a plateau in the velocity dispersion profile (for a cold 
collapse). Also, the apparent mass-to-light ratio for our clus- 
ters with DM is generally close to the case of a purely stellar 
cluster. (The special case of a hot collapse inside a DM halo 
which produces inflated values of T can be mistaken for a 
cluster being at its last stage of disruption by the tidal forces.) 

We argue that increasingly eccentric isodensity contours 
observed in the outskirts of some GCs could be created by 
a stellar cluster relaxing inside a triaxial DM minihalo, and 
not by external tidal fields as it is usually interpreted. In- 
deed, cosmological DM halos are known to have noticeably 
non-spherical shapes; a stellar cluster relaxing inside such a 
halo would have close to spherical distribution in its denser 
part where the stars dominate DM, and would exhibit isoden- 
sity contours of increasingly larger eccentricity in its outskirts 
where DM becomes the dominant mass component. 

It is also important to remember that few Galactic GCs 
show clear signs of a "tidal" cutoff in the outer density profiles 
(Trag er"et alJl995h . The "tidal" features of an opposite nature 
— "breaks" in the outer parts of the radial surface brightness 
profiles in some GCs — are often observed at or below the in- 
ferred level of contamination by foreground/background ob- 
jects, and could be in many cases an artifact of the background 
subtraction procedure which relies heavily on the assumption 
that the background objects are smoothly distributed across 
the field of view. A good example is tha t of D raco dwarf 
spheroidal galaxy. Ilrwin & Ha tzidimitriou ( 1995) used sim- 
ple non-filtered stellar counts from photographic plates fol- 



lowed by a background subtraction procedure to show that 
this galaxy appears to have a relatively small value of its ra- 
dial density "tidal cutoff" r, = 28.3 ± 2.4 arcmin and a sub- 
stantia l population of "extratidal" stars. lOdenkirchen et all 
(2001) used a more advanced approach of multi-color filter- 
ing of Draco stars from the Sloan Digital Sky Survey im- 
ag es and achieved much higher signal-to-noise ratio than 
in Ilrwin & Hatz idimitrioul i ll 9951) . New, higher quality re- 
sults were supposed to make the "extratidal" featu res of 
Draco much more visible. Instead, Odenkirchen et al. II2001 
demonstrated that the Draco's radial surface brightness profile 
is very regular down to a very low level (0.003 of the central 
surface brightness), and suggested a larger value for the King 
tidal radius of r, ~ 50 arcmin. 

We argue that the qualitative results presented in this pa- 
per are very general, and do not depend much on the the fact 
that we used MS04 model to set up the initial non-equilibrium 
stellar core configurations, and on the particular values of the 
model parameters (such as cr,.*, and x). As we discussed 
in Section[3] the appearance of "tidal" or "extratidal" features 
in our warm and cold collapse models is caused by two rea- 
sons: (1) at radii r > r m the potential is dominated by DM, 
whereas in the stellar core the potential is dominated by stars 
from the beginning till the end of simulations, and (2) the col- 
lapse is violent enough to eject a fraction of stars beyond the 
initial stellar cluster radi us. 

As we showed in § 12.41 for the initial stellar density 
p,> = 14 M Q pc" 3 (from MS04), the whole physically plau- 
sible range of stars-to-DM mass ratios x an d the GC forma- 
tion redshifts z satisfy the above first condition. For warm 
and cold collapse, this condition can be reexpressed as S = 
m*/mDM[r*] > 1- One can easily estimate S for other values 
of p i)it . 

The second condition is more difficult to quantify. In MS04 
we showed that collapsing homogeneous isothermal spheres 
produce extended halos for any values of the initial virial ra- 
tio v (except for v > 2.9 systems whic h are too hot to form 
a bound cluster in the absence of DM). iRov & Perezl ( 120041) 
simulated cold collapse for a wider spectrum of initial cluster 
configurations, including power-law density profile, clumpy 
and rotating clusters. In all their simulations an extended halo 
is formed after the initial violent relaxation. It appears that 
in many (probably most) stellar cluster configurations, which 
are not in detailed equilibrium initially, our second condition 
can be met. 

Of course, not every stellar cluster configuration will re- 
sult in a GC-like object, with a relatively large core and 
an extended halo, after the initial violent relaxation phase. 
Spitzer & Thuarj i 19721) demonstrated that a warm (v = 0.5) 
collapse of a homogeneous isothermal sphere produces clus- 
ters with large cores (their models D and G). Adiabatic col- 
lapse of homogeneous isothermal spheres produces clusters 
which surface density profiles are very close to those of dy- 
namically young Galactic GCs (MS04). Ir"ov & Perezl (120041) 
concluded that any cold stellar system which does not contain 
significant inhomoge neities relaxes to a la rge-core configu- 
ration. The results of Rov & Perez] ll2004l) can also be used 
to estimate the importance of adiabaticity for core formation. 
Indeed, their initially homogeneous models H and G span a 
large range of v, and include both adiabatic cases {v > 0.16 
for their number of particles N = 3 x 10 4 , from eg. 1161 ) and 
non-adiabatic ones (v < 0.16). It appears that all their mod- 
els (adiabatic and non-adiabatic) form a relatively large core. 
The issue is still open, but it appears that the adiabaticity re- 
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quirement (our eg. 1161 ) is not a very important one for our 
problem. (Though this requirement is met automatically for 
real GCs for the values of p, » and <t, * derived in MS04.) 

An important point to make is that the simulations pre- 
sented in this paper describe the collisionless phase of GC 
formation and evolution, and cannot be directly applied to 
GCs which have experienced significant secular evolution due 
to encounters between individual stars. In MS04 we showed 
that such collisionless simulations of purely stellar clusters 
describe very well the surface brightness profiles of dynam- 
ically young Galactic GCs (such as NGC 2419, NGC 5139, 
IC 4499, Arp 2, and Palomar 3 - see Fig. 1 in MS04). In 
Paper II we will address (among other things) the issue of 
long-term dynamic evolution of hybrid GCs. We will demon- 
strate that, at least for the warm-collapse case, secular evolu- 
tion does not change our qualitative results presented in this 
paper. 

In the light of the results presented in this paper and the 
above arguments, we argue that additional observational evi- 
dence is required to determine with any degree of confidence 

APPENDIX 

A. DISTRIBUTION FUNCTION FOR BURKERT HALOS 
The dimensionless potential * of Burkert halos with the density profile given by equations Q-(|8j is as follows: 



if GCs have any DM presently attached to them, or if they 
are purely stellar systems truncated by the tidal field of the 
host galaxy. A decisive evidence would be the presence of 
obvious tidal tails. A beau tiful example is given by Palomar 5 
( Odenkirchen et al. 2003), where tidal tails were observed to 
extend over 10° in the sky. 

Even if a GC cluster is proven not to have a significant 
amount of DM, it does not preclude it having been formed 
originally inside a DM miniha lo. I n "semi-consistent" simu- 
lations of Bromm & Clarke (2002) of a dwarf galaxy forma- 
tion, proto-GCs were observed to form inside DM minihalos, 
with the DM being lost during the violent relaxation accompa- 
nying the formation of the dwarf galaxy. Unfortunately, their 
simulations did not have enough resolution to clarify the fate 
of DM in proto-GCs. In Paper II we will address the issue of 
the fate of DM in hybrid proto-GCs experiencing severe tidal 
stripping in the potential of the host dwarf galaxy. 



<&=!• 



-1(1 +/T 1 )[arctan J R-ln(l +R)]+-(l -/T')ln(l +R 2 ) 

TT 2 



(Al) 



Here R = r/r s is the radial distance in scale radius units, and ^ is in TT 2 Gp^^r 2 units. The potential W is equal to 1 at the center 
of the halo, and zero in the infinity. 

Phase-space distribution function for a Burkert halo with an isotropic dispersion tensor can be derived through an Abel trans- 
form (Bin nev & Tremainel 19871 p. 237; we skip the second part of the integrant, which is equal to zero for Burkert halos): 



F(E): 



1 



d 2 p <i* 



(A2) 



v^Wo d* 2 (E-V) 1/2 ' 

Here E is the relative energy in the same units as and p is the density in po.i units, which results in (7rr s ) 3 p^G 3 / 2 units for F. 

The function d 2 p/d^ 2 for a Burkert halo has the same asymptotic behavior for R —* and R — > oo as the model I of Widrowl 
(2000), which has the density profile p oc (1 + /?)~ 3 : 



d 2 p 



(R -> 0) oc R~ 



d 2 p 



(R — * oo) oc 



d^ 2K " ' c/* 2 (-In*) 3 ' 

This allowed us to use the following analytic fitting formula of Widrow (2000) for the distribution function of Burkert halo: 



(A3) 



F(E) = F E y2 ( l-ET 1 



-\nE 
l-E 



Here P = J^iPi^' i s a polynomial introduced to improve the fit. 

We solve equation ( IA2t numeri cally for the interval of radial distances R = 10~ 6 . 
the fitting coefficients in equation (IA4> : 



(A4) 



10 6 . We obtained the following values for 



F = 5.93859 x 10" 4 ; ^ = -2.58496; p x =-0.875182; 

p 2 = 24.4945; p 3 =-147.871; p A = 460.351; (A5) 

p 5 = -747.347; p b = 605.212; p n =-193.621. 

We checked the accuracy of the fitting formula (IA4> with the above values of the fitting coefficients by compa ring the an alytic 
densi ty profile of the Burkert halo (eq. Q) with the density profile derived from the distribution function ( Bi nnev & Tremainel 
CHI p. 236): 



p(R) = iVItt 2 I F(E)(^-E) l/2 dE. 
Jo 



(A6) 



The deviation of p(R) in equation iA6\ . with F(E) given by equations ( IA4> - (IA5> . from the analytic density profile was found to 
be less than 0.7% for the interval of the radial distances R = 10~ 6 . . . 10 6 . 
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FIG. B 1 . — Scheme explaining the projection procedure. The stellar cluster is marked with its center O and thin line circular boundary. The axis Z is arbitrarily 
chosen to point toward observer. We show a sphere containing the i-th particle Pj, located at the distance r, from the center of the cluster, as a thick line circle. 
The range of projected distances Ri ...i?2 is shown as two vertical shadowed bars. In 3D, it corresponds to space between two concentric cylinders of radii R[ 
and R2 with the axis OZ as their central axis. These two cylinders cut segments out of the sphere containing the particle Pj, which are marked as four regions A, 
B, C, and D. In 3D, these four regions correspond to two segments of the sphere. For an arbitrary rotation of the cluster around its center O, the probability to 
find the particle P, within the bin R\ . . .R2 is equal to the ratio of the surface area of the segments to the total surface area of the sphere 4nrj. 



B. PROJECTION METHOD 

For iV-body models of spherically symmetric stellar systems, such as GC models presented in the paper, one can produce radial 
profiles of projected observable quantities (such as surface mass density or velocity dispersion) by averaging these quantities over 
all possible rotations of the cluster relative to the observer. Radial profiles created in this way preserve maximum information, 
which is very important for studies of the outskirts of the clusters and for clusters simulated with relatively small number of 
particles. 

Let us consi der a n i-th stellar particle P,- in a spherically symmetric cluster, located at the distance r, from the center of the 
cluster (Figure IbTI . We want to find its averaged contribution to one projection bin, corresponding to the range of projected 
distan ces R i . . .R 2 . For an arbitrarily rotated cluster, the probability w, to detect our particle inside the bin (the shaded areas in 
Figure lBlTi is equal to the fraction of the surface area of the sphere, containing the particle f,, cut out by two concentric cylinders 
with the radii R\ and R2- It is easy to show that this probability is 

TO, n<Ri; 
Wi = J [l-0rVr ; ) 2 ] I/2 , Ri^n<R 2 ; (Bl) 

[[l-(^/r,) 2 ] 1/2 -[l-(/? 2 A,.) 2 ] 1/2 , r t >R 2 . 

The projected surface mass density ( for the bin R\...R 2 averaged over all possible rotations of the cluster is 

1 N 

C(RuR 2 )= = t-ymiWi, (B2) 

S ' 7T(fl 2 -/? 2 )^ ' 



where N is the total number of stellar particles in the cluster and m, is the mass of the i-th particle. The averaged value of a 
projected scalar quantity Q is 

N I N \ _1 

Q(RuR2) = J^w i Q i (l><j . ( B3 > 

so the probability w,- serves as a weight for scalar variables. For example, to calculate the square of the projected three-dimensional 
velo city dispersion ct^d, corresponding to the projected distances range R\ . . .R2, one has to set Q, = (V xi + V?,- + V. 2 ;) in equa- 
tion jB3> , where (Y x j,V y j,V Z i) are the velocity components of the i-th particle. Observationally, one can determine a^ D by mea- 
suring both line-of-sight velocities and two proper motion components for stars located within the interval R\ . . .R 2 of projected 
radial distances. 

One can show that the above method is superior to a straightforward projection onto one plane or three orthogonal planes by 
considering a contribution made by one particle to the bin R\ .. .R 2 . In case of the straightforward projection, there is a high 
probability that the particle will make no contribution to the projection bin. In our method, every particle with r, > Ri will make 
a contribution to the bin, which will make the radial profile of the proj e cted quantity significantly less noisy. Central value of the 
projected quantity is obtained by setting Ri to zero in equations (IB1I — dB3t . In this case, all particles make contribution to the 
averaged value of the quantity. 

The situation with vector quantities is a bit more complicated. We designed a simplified procedure which is based on the 
assumption that the bin sizes are small (so the number of bins is large). Each particle with an associated vector quantity (say, 
velocity vector) is first rotated around the OZ axis to bring the particle in the plane XOZ. Next, the particle is rotated around the 
axis OY to make the projected radial distance of the particle equal to the radial distance of the center of the bin R = {R\+R2)/2. 
(The observer is assumed to be at the end of the axis OZ.) 



14 



Mashchenko and Sills 



The detailed procedure is as follows. For a particle with coordinates (x,y,z), after the two rotations the associated vector 
(V x ,V y ,V z ) is transformed into vector (V" ,V" ,V") with the following components: 

V"=V^sinip\ V" = V fV shV, V" =V^cosip' . (B4) 
The parameters V' K „ V v , ip', and cp' can be obtained after a series of the following calculations: 

r = (x 2 +y 2 + z 2 ) l/2 ; r xy = {x 2 +y 2 ) 1 l 2 ; V xy = (V 2 + V 2 ) 1 ^ 2 ; 

sincp = V y /V xy ; coscp = V x /V xy ; sina = y/r xy ; 

ip' = ip— a; 

V^ = (V 2 ,cos 2 (p' + V 2 )^ 2 ; 

sin ip = V xy cos p'/V^; cos tjj = V Z /V X7 ; sin 9 = r„/ r; 
t/j' = ip - 9 + arcsin x. 

Here x is equal to R/r if R/r < 1 and is equal to 1 otherwise; r, r xy , V xy , cp, a, 9, and ip are intermediate variables. 

If vector (V x ,V y ,V z ) is velocity vector, then substituting V" 2 in place of Qj in equation JB3I will give us the value of the square 
of the line-of-sight velocity dispersion a r , averaged over all possible rotations of the cluster, for the bin R[ . . ./?2- Similarly, 
replacing Qj with V v " 2 and V'.' 2 will produce the averaged values of the square of the radial and tangential components of the 
proper motion dispersion, respectively. 



cosa = x/r xy ; 
cos 9 = z/r; 



(B5) 
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